clear; clc; close all;
% Lyapunov 指数 vs. z0

mu = 3.5;
k  = 0.16;
x0 = 0.1;
y0 = 0;

z0_range = linspace(-10, 8, 1000);

% 迭代步数设置
N_transient = 1000;   % 消除瞬态的迭代步数
nLE_steps   = 500;    % 用于计算 Lyapunov 指数的迭代步数

% 发散判定阈值
threshold = 1e6;

% 预分配存储 LE 值的数组，每行对应一个 z0 值，3列分别为 LE1, LE2, LE3
LEs_values = nan(length(z0_range), 3);

% 对每个 z0 值进行计算
for i = 1:length(z0_range)
    z0 = z0_range(i);
    state = [x0, y0, z0];
    diverged = false;
    % 瞬态迭代
    for t = 1:N_transient
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
        if any(abs(state) > threshold)
            diverged = true;
            break;
        end
    end
    
    if diverged
        % 发散的情况，赋值为 NaN
        LEs_values(i, :) = [NaN, NaN, NaN];
    else
        % 计算 Lyapunov 指数，使用你的 LEs.m 函数
        LEs_values(i, :) = LEs(state, mu, k, nLE_steps);
    end
end

figure;
plot(z0_range, LEs_values(:,1), 'r', 'LineWidth', 1.5); hold on;
plot(z0_range, LEs_values(:,2), 'g', 'LineWidth', 1.5);
plot(z0_range, LEs_values(:,3), 'b', 'LineWidth', 1.5);
xlabel('$z_0$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('LEs', 'Interpreter', 'latex', 'FontSize', 14);
title('Lyapunov 指数谱');
legend('LE1','LE2','LE3','Location','Best');
grid on;
